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Abstract 

We consider a generalisation of Ulam's method for approximating invariant den- 
^ sities of one-dimensional chaotic maps. Rather than use piecewise constant polyno- 

^ f-H mials to approximate the density, we use polynomials of degree n which are defined 

^ by the requirement that they preserve the measure on n + 1 neighbouring subinter- 

^ vals. Over the whole interval, this results in a discontinuous piecewise polynomial 

I— I approximation to the density. We prove error results where this approach is used to 

^ approximate smooth densities. We also consider the computation of the Lyapunov 

^ exponent using the polynomial density and show that the order of convergence is 

r one order better than for the density itself. Together with using cubic polynomi- 

als in the density approximation, this yields a very efficient method for computing 
i/^ highly accurate estimates of the Lyapunov exponent. We illustrate the theoretical 

^ 5 findings with some examples. 

Tl 1 Introduction 

> 

•i-H 

^ Suppose that a one-dimensional piecewise smooth map g \ I ^ I on compact interval 

^ / C M has a unique chaotic attractor with an absolutely continuous invariant measure 

/i. We consider piecewise polynomial approximations to the density d associated with /i. 
To this end, we define a partition of the interval / by dividing it into a number of equal 
subintervals Ii = i — 1,2, . . . , A/', each of length h — {x^ — xq)/N. For each 

subinterval, we define 

rrii — I dji^ i = 1, 2, . . . , A^. 
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Thus, rrii is the mass of the interval li with respect to the invariant measure /x. Clearly 
we have the property that 

N 
i=l 

since /i is a probability measure. We want to approximate the density d oi ji given that 
the only information that we have about /i is the masses m^, i = l,2,...,A^. 

For approximating invariant densities, Ulam's classical method [39] consists in using 
piecewise constant approximations on each subinterval with the constants being found 
by computing a fixed point (i.e. an eigenvector at the eigenvalue 1) of a (discretized) 
Markov operator, the so called Frobenius- Perron operator or more generally transfer op- 
erator. This operator describes how measures (and densities) are pushed forward under 
the dynamics of the map g. Ulam's method has been proved to be convergent for piece- 
wise expanding maps [34]. Since then much effort has gone into analysing the method 
[2, 7, 12, 14, 15, 20], computational questions [13, 27], extending the proof to more general 
classes of transformations [16, 22, 23, 24, 28, 29, 30, 35] or extending the method to higher 
dimensional systems [8, 9, 10, 18, 21, 31, 32]. 

For non-smooth invariant densities, a piecewise constant approximation may be a 
good choice. However, in the smooth case (for example if the map is perturbed by white 
noise), it would be better to use a higher order approximation. Higher order polynomial 
approximations to the invariant density have been obtained by Ding and coworkers [11, 15, 
17, 19]. They used a Galerkin method for finding a piecewise polynomial approximation 
to the density but did not go on to consider the computation of the Lyapunov exponent. 

Increasingly, numerical methods are being developed which incorporate and preserve 
features of the problem being considered. One well-known example of this is geometric 
integrators for ODE's which encode properties of the ODE's, such as a symplectic struc- 
ture or conservation of an invariant, into the numerical method [26]. Our approach to the 
problem of computing an invariant density is based on this philosophy. The essence of the 
density is that the measure should be preserved on any subinterval and so we construct 
a discontinuous piecewise polynomial basis which is defined locally by the requirement 
that the measure is preserved on neighbouring intervals. Combining many such polyno- 
mials over the whole interval, the coefficients of the approximate invariant density with 
respect to this basis can again be found by computing a fixed point of the associated 
discrete Frobenius-Perron operator. Using this new basis, which turns out to lead to a 
Petrov- Galerkin discretisation of the Frobenius-Perron operator, we improve upon the 
convergence rates in [11, 15] for piecewise quadratic approximations. In addition, we 
consider the computation of the Lyapunov exponent using the invariant density and give 
corresponding error estimates. It turns out (see Theorem 5.1) that using our measure pre- 
serving basis, one actually gains one order in the convergence of the Lyapunov exponent. 
Using cubic polynomials, this leads to a highly efficient method for computing estimates 
of the Lyapunov exponent. 

An outline of the paper is as follows: In Section 2, we derive the measure-preserving 
polynomial basis and perform an error analysis for a density that is approximated by such 
a polynomial. Section 3 considers the problem of integration when the density is replaced 
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by a polynomial while Section 4 gives an error analysis when the density is found as a 
solution of the discretised Frobenius-Perron fixed point equation. The computation of 
the Lyapunov exponent using the polynomial approximation to the density is considered 
in Section 5, and the errors in the Lyapunov exponent are derived. These results are 
illustrated with some examples. Finally, in Section 6, two extensions of this work are 
described. 



2 Measure-Preserving Polynomial Approximations to 
the Density 

Our aim is to construct a piecewise polynomial approximation to an invariant density. To 
do this, we start by considering the problem of obtaining a polynomial approximation to 
the density function over a subset of n + 1 adjacent intervals [x^, Xi^n+i] = ^]=o^i+j+i for 
some i. As with Gaussian quadrature, we prefer to work with the standard interval [—1,1] 
and so we perform a change of variables to map the interval [x^, Xi^n+i] onto [—1,1] using 

The corresponding inverse transformation is given by 

x{t) = -{n + l)h{t + l) + Xi. (2.2) 

We transform the points x^+j onto the interval [—1,1] by 

t, = t{x,^j) = -1 + J = 0, 1, . . . , n + 1, (2.3) 

assuming that the points are evenly spaced. We also define a new density function D{t) 
which satisfies 

D{t) dt = d{x) rfx, 

which implies that 

dr 1 

D{t) = -d{x{t)) = -{n + l)hd{x{t)), (2.4) 
using (2.2). Using this definition, we find that 

/ f{x)d{x) dx= I F{t)D{t) dt, (2.5) 

J Xi J —1 

where F{t) = f{x{t)). Defining 

Mj = rrii^j^i, j = 0, . . . , n, (2.6) 

we obtain 

Mj m^+j+i / d{x) dx = D{t) dt^ j = 0, . . . , n. 
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We now use the n + 1 values Mj, j = 0, . . . , n to obtain a polynomial Pn{t) G 11^ which 
approximates the density function D{t) which can be used in the integration formula 
(2.5) with a higher order Gaussian quadrature method to obtain higher accuracy for the 
integral. The criterion which we use to define the polynomial is that it should preserve 
the measure on each subinterval, that is 

r^\n{t) dt = Mj, J = 0,l,...,n. (2.7) 

We will consider two diflFerent methods for determining the measure-preserving polynomial 
Pn{t)^ but we first establish uniqueness. 

Lemma 2.1. The polynomial Pn{t) G 11^ satisfying (2.7) is unique. 

Proof. Assume that Pn{t) and qn{t) are both polynomials that satisfy the conditions (2.7) 
and define 

hn{t) =Pn{t) - Quit). 

Then hn{t) G Un and satisfies the conditions 

r^\n{t) dt = 0, J =0,l,...,n. 
Jtj 

These conditions imply that, for each value of j, there is at least one point tj G (tj,tj+i) 
such that hn{Tj) = 0. Thus, in total, there are at least n + 1 distinct points in the interval 
(—1,1) at which hn vanishes. The only polynomial of degree n that has at least n+1 
distinct zeros is the zero polynomial, and so hn{t) = 0, giving Pn{t) = qn{t). □ 

There are two approaches which can be used to determine the polynomial Pn{t)^ both 
of which will be useful later. 

2.1 Constructing the measure-preserving polynomial by inter- 
polating the measure 

We define the measure u associated with the density function D{t) by 

V{t)= j\{T)dT, (2.8) 

and we note that 

i-i 

v{tj) = 5^ M,, i = 1, 2, . . . , n + 1, z/(to) = 0. (2.9) 

Strictly speaking, v is not a probability measure since J^-^ rfz/ 7^ 1, but it is a part of the 
probability measure ji. From (2.9), we have n + 2 values of the function iy{t) which we 
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can interpolate using a polynomial G H^+i. The derivative of this polynomial will 

give an approximation to the density function which satisfies (2.7). Thus, we define 

where prime denotes differentiation with respect to t. Using the Lagrange form of the 
interpolating polynomial [5] , we thus have that 

j=l 3=1 k=0 k=0 j=k+l 

where 



n . 



Note that there is no term in the definition of with j = since z^(to) = 0. Thus, 

one representation of the polynomial which satisfies the conditions (2.7) is 

n 72+1 
k=0 j=/c+l 

2.2 Constructing the measure-preserving polynomial using an 
appropriate basis 

The derivatives of the basis functions Lnj{t) in the previous method are clearly quite 
messy to evaluate. An alternative approach is to seek the polynomial Pn{t) in the form 

n 
k=0 

The basis functions £n,k{t) can be found by substituting this form of the polynomial 
into the conditions (2.7) and equating coefficients of M^- This implies that for each 
= 0, 1, . . . , n, we have the n + 1 conditions 

iuAt) dt={ ^ J = 0,l,...,n. (2.11) 

These equations represent a system of linear equations for the n + 1 coefficients in the 
polynomial in,k{t). 

Of course this is only a different representation of the same polynomial found using 
the measure and so 

j=k+l 

The basis functions defined by the conditions (2.11) for some low values of n are listed in 
Table 1. The three basis functions for the case n = 2 are shown in Fig. 1. Some properties 
of these basis functions can be determined. 
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n = 


ko{t) = I 


n = 1 


=^i,o(-^) 


n = 2 


£2,o(t) = 1^(27^2 -18t-l) 
^2,1 W = |(-27t2 + 13) 

4,2(t) = 4,0(-t) 


n = 3 


^3,oW = |(-16t^ + 12^2 + 2t - 1) 
l^ -^it) = 1(48^3 - 12^2 - 30t + 7) 

4,2(t)=4,i(-^) 



Table 1: Basis functions for the measure-preserving polynomial approximation to the 
density. 

Lemma 2.2. The basis functions £n,k{t) have the following properties: 

(i) in,n-k{t) = 

(ii) J JnAi) dt = I; 

(iii) X^^«.'^(^) = ^' 

k=0 

(iv) The function £n,k{t) has n simple roots with one root in each of the intervals 
(tj^tj^i) for j ^ k. 

Proof. Results (i) and (ii) follow directly from the definition of the basis functions. For the 
third result, consider the case that D{t) = c for some constant c. Then Mi = 2c/ {n + 1), 
i = 0, . . . , n and Pn{t) = D{t) = c. Substituting these into (2.10) gives the required result. 
For the fourth result, note that the condition that the integral of in,k{t) is zero over the 
subinterval [tj^tj^i] implies that there must be at least one point in (tj, tj+i) at which the 
function is zero by the first Mean Value Theorem for Integrals. Since one root in each 
interval gives us the maximum possible number of roots, then there must be precisely one 
root in each interval. □ 
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-2.0^ 



Figure 1: The three basis functions for n = 2. 



2.3 Error analysis 

Having determined methods for finding a polynomial approximation Pn{t) to the density 
function D{t)^ we must now consider the errors in this approximation. Thus, we write 



(2.12) 



for some error function en{t). 

Using the first approach of interpolating the measure, the error in the interpolating 
polynomial is given by 



where 



n+l 



it) 



(n + 2)! 11^^ 



(2.13) 



for some /^(t) G (—1,1), assuming that iy{t) G 1,1] (see [5]). However, to obtain 

the error in the density requires differentiating 6:^+1 (t) which does not lead to a nice form 
for en{t). Thus, we also derive the error based on the alternative approach of determining 
the polynomial in terms of the basis functions in,k{t)- 

Integrating (2.12) over one subinterval and using (2.7) gives that 



en{t) dt = 0, j = 0, 1, . . . ,n. 



(2.14) 



By the first Mean Value Theorem for Integrals, this implies that there is at least one 
point Tj G (tj, tj+i) such that en{Tj) = 0, j = 0, 1, . . . , n and so there are at least n + l 
distinct points in the interval (—1,1) at which the polynomial interpolates the density. If 
we choose precisely one such point in each subinterval, then this interpolation problem 
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defines a unique polynomial in 11^ which must be the polynomial Pn{t)' Thus, regarding 
the problem as an interpolation problem, we can express the error function as 

(n + l)! n*'""^'- 

where ^{t) G (—1, 1) assuming that D{t) G 1, 1] (see [5]). Note that the points tj 

depend on the particular function D{t). 

Since there is a simple relationship (2.8) between iy{t) and D{t)^ there must also be a 
similar relationship between the errors which is given by 

en+i{t) = [ en{T)dT. (2.16) 



3 Integration 

There are two diflFerent types of integral that will be needed later which involve the density 
and so we consider them both at this stage, using the polynomial approximation to the 
density, and also provide an error analysis. 



3.1 Case I 

When computing the density using the Frobenius-Perron equation, integrals of the form 

d{x) dx = D{t) dt 

will be required, where = t(x^), tr = t{xr) and [t^,tr] ^ ["l^l]- We assume in this 
section that the measures rui on the intervals i = 1, . . . , are known exactly and we 
consider only the errors associated with the piecewise polynomial approximation to the 
measure. Using the polynomial approximation to the density then gives the following 
result. 

Theorem 3.1. If D{t) G then 



ti 

tr 



D{t) dt = I pn{t) dt + 0(/i^+2), (3.1 



r \D{t)-pn{t)\ dt = 0{h^^'). (3.2) 

Jte 

Proof. Integrating the polynomial approximation with error given by (2.12), we obtain 



ntr rtr rtr 

I D{t)dt= I pn{t)dt+ / Bnit) dt. 

J t£ J t£ J t£ 
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Thus, for (3.1), we need to estimate 

rtr 

Similarly, for (3.2), we need to estimate 



[\n{t) dt. (3.3) 



^ \en{t)\ dt. (3.4) 

From (2.15), we see that the function en{t) involves D^^^^\^{t)). Now from (2.4), we 
know that 

D{t) = chd{x{t)), 
where c = |(n + 1). Since dx/dt = ch^ we therefore note that 

i^("+^)(t) = (c/i)"+2rf("+^)(x(t)), 

and so 

i^("+i)(e(t)) = 0(/i"+2). 
Thus, the integrals (3.3) and (3.4) are both 0(/i^+^) as required. □ 

Using (3.1), we now obtain 

u 



PXr rtr 

/ d{x) dx = D{t) dt 



\n{t) dt + Oih""^^) 

tr ^ 



k=0 

= 5j^m,+,+i r4,^(t) rft + 0(/i^+2), 

k=0 

where we used (2.6) in the final step. Clearly, in practice we simply drop the error term 
when computing this integral. 

3.2 Case II 

The second type of integral that we will need to evaluate is given by (2.5). As in the 
previous case, we again assume that the measures rrii on the intervals i = 1, . . . , are 
known exactly. Then we have that 

J F{t)D{t)dt = J F{t)pn{t)dt + J F{t)en{t) dt. (3.5) 
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We first consider tlie error term in (3.5). For later use, we define tlie function 

n+l 

<j)n+2{t) = l[{t-tj), 
j=0 

where the evenly spaced points tj are given by (2.3). 

Theorem 3.2. If n is odd, F G C^[— 1, 1] and D G C"'"'"^[— 1, 1] then the error term in 
(3.5) is 



F{t)en{t) dt 



F'\a)D^^+^\5) ^ F'{a)D^^+^^{-f) 



(n + 2)! 



(n + 3)! 



j^n+2{t)dt, (3.6) 



where a,S,jE (—1, 1). 

If n is even, F e C^[— 1, 1] and D e C""'""^[— 1, 1] then the error term in (3.5) is 



F{t)en{t) dt = 



F"{a)D^^+^\5) ^ F'(«)D("+2)(7) 



(n + 2)! 



(a - 1) / (t>n+2{t) dt 



(n + 3)! 



F(OZ}("+^) (77) 
(n + 2)! 



4>n+2it) dt, 



(3.7) 



where a, 5, 7, G (—1, 1). 



Proof. We analyse the error by integrating by parts and using the error in the approxi- 
mation to the measure v{t) given by (2.13). Integrating by parts and using (2.16) gives 



F{t)en{t) dt = [F{t)en+i{t)]\ - / F'(t)5.+i(t) dt 



The boundary terms vanish since 6:^+1 (—1) = by (2.16) and 6:^+1 (1) = by (2.16) and 
(2.14). Thus, 



1 



F\ty^+^\K{t))<P^+2{t) dt. 



(3. 



This error term is very similar to that for simple Newton-Cotes integration formulae [37] 
since the tj's are evenly spaced points, except for the extra term F'{t) which complicates 
the analysis. 

We first consider the case when n is odd. In this case, 0^+2 (^) is an odd function since 
the points tj are evenly spaced, and this results in an error which is of higher order than 
would be expected, as in the case of Newton-Cotes formulae. 

We define 



V) dr. 



1 
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Since 0^+2 (^) is an odd function, it follows that ipn+sit) is of one sign [37], and that 

^Pn+3{l) = 0. (3.9) 
Integrating (3.8) by parts, the boundary terms again disappear using (3.9) and so 

F'it)u^-^'\Kim^+2it) dt = - £^ ^ {F'{t)u^^+'\K{t))) ^„+3(t) dt. 
Since '0„+3(t) is of one sign, the Mean Value Theorem for Integrals can be used, giving 



d 



F'{t)u(-+'\K{t))<t>^^,{t) dt = -- {F'{ty-^'\K{t)))l^^ / >n+3(t) dt 



dt 



{F'{t)u(^+'HMmLa / t<j>^+2{t) dt,{3.10) 



where a G (—1, 1) and the second step involved another integration by parts. 
It can be shown [36] that 



d 



,(n+2) , 



,(n+3) 



m)) 



(3.11) 



dt \ \ ' ' n + 3 

for some /3{t) G (—1, 1). Thus, substituting (3.10) into (3.8) and using (3.11), we obtain 



F{t)en{t) dt = - 



F"{ay^+^\5) ^ F'(«)z^("+')(7) 



(n + 2)! 



(n + 3)! 



/: 



t(l)n+2{t) dt, 



where 5 = tv{a) and 7 = /3{a). Now iy\t) = D{t) and so substituting for u gives (3.6) as 
claimed. 

When n is even, the proof is more complicated, as it is for the standard Newton-Cotes 
result. We first break the integral in (3.8) into two, giving 

F\ty-^'\^{t))cl^r^^,{t) M = y^%^(t)z.(-+2)(/.(t))0.+2(t) 

+ rFXt)z.(-+2)(/.(t))0.+2(t)rft. (3.12) 

For the second integral, 0^+2 (i) is of one sign on the interval [t^, 1] and so the Mean Value 
Theorem for Integrals gives that 

r F'{t)iy^^+'\K{t))<t>„+2{t) dt = F'{ey-+'\K{e)) f 4>^+2{t) dt, (3.13) 

for some 9 G (t^i, !)• 

We rewrite the first integral in (3.12) as 



F'{ty^+'\K{t))<t>^+2{t) dt 



F\ty^+'\K{t)){t - l)(t>^+,{t) dt, 
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and we note that 0^+1 (t) is an odd function about the midpoint of the interval [— 
Thus, using the result (3.10) for n odd then gives 

F'{ty-+'\^m^^,{t) dt = {F'{ty-^'\Kmt - i)}L_ x 



t<j)n+l{t) dt 



d 



-{F'{t)u^-^-\^{t))]l^Ja-l) 



where a G (— Now since (f)n+iit) is odd about the midpoint of the interval, 



t<^„+l(t) dt= {t- l)0n+i(t) dt = / 4>n+2{t) dt. 



(3.15) 



Combining the integrals (3.13) and (3.14), and using (3.15), we therefore obtain 



F'{t)u^-^'\K{t))<i>r.+2{t) dt 



j^{F'{t)u^-^^\K{t))}l^Ja-l) 



+ F'(a)i/("+2)(fi:(a)) 

+F\9)u^'^^'\k{9)) [\r.+2it)dt. 

Now the two integrals have the same sign [38] and so 

f_F'{t)v^-^'\Kmn^,{t)dt = f^{F'{t)u^-+'\K{t))}l^^x 

rtn 

(« - 1) / (l)n+2{t) dt 



(t) dt 



(t) dt+ / (j)n+2{t) 
-1 Jt„ 

{F'{ty-^'\K{t))]l^^x 



dt 



dt 



a -I) I (t)n+2{t) dt 



+F\0i^(-^'\<0)l\+2{t)dt, 
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for some G (— 1, 1) (and satisfying a < ^ < 9). Expanding the derivative term and using 
(3.11) gives 



F'{t)iy'^^+'\K{t))<P^^2{t) dt 



F'(a)z/("+3)(/3(a)) 



n + 3 



tn 



x(a - 1) / (j)n+2{t) dt 



+F\0'^^''^'\<0)j\+2{t)dt, 



where /3{a) G (—1, 1). Substituting back into (3.8) then gives 



F{t)en{t) dt 



F"(a)z^("+2) («;(«)] ^ F'(a)i/("+3)(^(«)) 



(n + 2)! 



(n + 3)! 



X 



ia-1] 



(j)n+2{t) dt 



(n + 2)! 



4>n+2{t) dt, 



Finally, substituting for p using u'{t) = D(t) gives the stated result, with 5 = /«(«), 
7 = /5(q;) and ?] = /«(^). □ 

These results can be extended to the case of evaluating an integral over the whole 
interval /. 

Theorem 3.3. Let p%{x) he the discontinuous piecewise polynomial approximation to 
the density, assuming that the measures nii on the intervals li, i = 1, . . . , A^" are known 
exactly. 

If n is odd, f G C^i^I), d G C"'"'"^(/) and f(x) is not a constant, then 



f{x)d{x) dx = I /(a;)j9^(a;) dx + 0(/i"+2). 



If n is even, f G C^{I), d G C"~^^{I) and f(x) is not a constant, then 



f{x)d{x) dx = / f{x)p%{x) dx + 0(/i"+^) 



(3.16) 



(3.17) 



Proof. We first note that if f{x) is a constant, then the error term in both cases is zero 
by construction, and hence this case has been excluded. The results of Theorem 3.2 can 
be converted back to the original x coordinates in order to obtain error terms as a power 
of h. As in the proof of Theorem 3.1, we note that D{t) = chd{x{t)) and that dx/dt = ch^ 
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where c = |(n + 1). Using these results, we obtain, for n odd, 

Xi-\-n-\-l /*! 

f{x)d{x) dx = F{t)D{t) dt 

F{t)pn{t) dt+ [ F{t)en{t) dt 
1 J-1 
1 

F{t)pn{t) dt + 0{N'^^). 

1 

Similarly, for n even, we have 

rxij^n+i rl 

/ f{x)d{x) dx= F{t)pn{t) dt + 0(/i^+^). 

Jxi J-1 

We note that in the first term of (3.7), the term a — 1 can be expressed, using (2.1), as 

^ 1 o , 2(a-x^) 2{a-Xi+n+i) n(Mh\ 

01 — ^ = —2 + = — 7 = Oil n), 

{n + l)h {n + l)h ^ ' ^ 

where a = x{a) G (x^,x^+^+i), and so both terms in (3.7) are 0{h^^^). 

Integrating over the whole interval / requires the sum of N/ (n + 1) such integrals. As 
usual, the error term for such a composite quadrature rule is one power of h less than 
for the simple rule as it involves the sum of the errors over the N /{n + 1) intervals, and 

= 0{h~^). This gives the stated results. □ 

Theorem 3.4. The results of Theorem 3.3 hold if f G C^{I) hut is piecewise (n odd) 
or piecewise (n even). 

Proof. Clearly for intervals li which do not contain a discontinuity in the derivative of 
/(x). Theorem 3.2 still holds. Thus, we need consider only those intervals where there 
is point of discontinuity of the derivative. In particular, we will consider only the case 
where there is a single point of discontinuity in any particular interval, but the results 
can easily be generalised to the case of multiple points of discontinuity. Thus, we assume 
that on some interval the function F{t) = f{x{t)) is given by 



m = 



Fi{t), -l<t<t* 
F2(t), t* <t<l 



where Fi{t*) = F2{t*), but F[{t*) ^ F'^(^*). The result obtained from the first integration 
by parts in the proof of Theorem 3.2 stiU holds since F is continuous. However, the 
resulting integral must be split into two and so we have 



^ 1 



II*' 

F(t)e„(t) dt 



(n + 2)! 



F;(t)z/(-+2)(«:(t))0„+2(t) dt^ 

1 

[' Fiit)u^^^'\K{t))4>n+2{t) dt 
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Ti 


Flrror T'prrn 




1 

2 
3 
4 


0{h^) 



Table 2: The order of the error term in (3.16) and (3.17) for different values of n. 

For n odd, the next step in the proof of Theorem 3.2 was to integrate by parts again, with 
the boundary terms vanishing. However, in this case, the boundary terms do not vanish, 
but integration by parts gives 

nt)en{t) dt = {F[{f) - F^{f)) v^-^'\n{f))Mn 

£f^{F^{t)i^^-^'\Kit)))^^^sit)dt . 

Now ipn+sit) is of one sign [37] and so, by the Mean Value Theorem for Integrals, the 
derivative term can be taken outside of each integral, giving two terms each of which 
are similar to that obtained in Theorem 3.2. Converting back to x coordinates, as in 
the proof of Theorem 3.3, we again find that the integral terms are 0{h^^^). However, 
the boundary terms, which were not present previously, are 0{h^^^). Now when we sum 
over all the intervals, we sum the integrals and this sum gives an error of 0(/i^^^) since 

= 0(/i~^), as previously. In the case of the boundary terms, we note that only intervals 
containing a discontinuity of the derivative of / have such boundary terms, and this is 
a fixed number which does not depend on h. Thus, the contribution of these boundary 
terms is also 0{h^^^) and so the total error is 0{h^^^) as previously. 

A similar analysis applies for n even. □ 

A summary of the order of errors for different values of n is given in Table 2, which 
shows a similar pattern to the error in Newton-Cotes integration. 

3.3 Gaussian quadrature 

In Case I, the integral on the right hand side of (3.1) is that of a polynomial, and so can be 
evaluated exactly. However, in Case II and for a general function F(t), the first integral 
on the right hand side of (3.5) may have to be evaluated numerically using Gaussian 
quadrature. When integrating a function G{t) on the interval [—1,1] using an m point 
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Gaussian quadrature rule, the error term is given by 



(2m)! 




where Pm{t) is the m^^ Legendre polynomial and rj G (—1,1) [5]. Converting back to 
the original x coordinates, we note again that each derivative of a function introduces 
a power of /i, since dx/dt = |(n + l)/i. The functions we integrate always involve the 
density function D{t) and the change of function from D{t) (which is approximated by 
Pn{t)) back to d{x) introduces another power of h. Thus, the error associated with a 
simple m point Gaussian quadrature rule is 0{h'^^^^). 

In Case II, as shown in the proof of Theorem 3.3, the approximation error for a single 
integral for n odd is 0{h'^^^). Thus, the error from the quadrature will be the same order 
as the approximation error if m = (n + 3)/2. Similarly for n even, the approximation 
error for a single integral is 0(/i^^^) and so the quadrature error will be the same order 
if m = (n + 2)/2. Thus, for n = 2/c — 1 or n = 2/c, the number of integration points used 
should be at least m = k + 1 to ensure that the error in the integration does not dominate 
the approximation error. 



4 Computing the Polynomial Density from the 
Frobenius-Perron Equation 

4.1 Discretisation of the Probenius-Perron equation 

Using the approach described in the previous section, we have reduced the description 
of the density associated with an iteration function g : I ^ I to Sl finite dimensional 
approximation which is characterised by the A^-dimensional vector m whose components 
are 

rrii = rf/i, i = 1, 2, . . . , A^. 



To determine this vector, we must solve the Frobenius-Perron equation given by 



dfi = dfi 

for any interval J C I. To obtain a set of determining equations for the vector m, we use 
the Frobenius-Perron equation with J = /^, i = l,2,...,A^ giving 



djji = I dfi^ z = 1, 2, . . . , A^. 



Now by definition, d/a = rrii. Similarly, the right hand side can be written as 

/ rf/x = / d{x) rfx. 
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where K is the number of preimages of the interval and for each k = 1^ . . . ^ K ^ g{xL^k) = 
Xi-i^ g{xR^k) = Xi. Now each of the preimages of the interval li may be contained within 
a single interval Ij for some j or it may contain some whole intervals and only parts of 
others. The evaluation of this type of integral was considered in Case I in the previous 
section. 

We note that when n = (piecewise constant polynomial approximations), this is the 
standard Ulam method for approximating the invariant density [39] . 



4.2 Convergence 

We now consider the convergence of this method as ^ oc. We define the approximation 
space 

= {/ : [0, 1] ^ [0, 1] : G n„, ^ = 0, n + 1, 2(n + 1), . . . , iV - (n + 1)} , 

of functions which are piecewise polynomials of degree n on n + 1 adjacent intervals from 
the given partition. Note that dim(A^) = and functions in are not necessarily 
continuous. 

In Section 2, for a given density rf, we constructed a function G which satisfies 



p%{x) dx = rrii = / d{x) dx^ i = 1, . . . , A^. 



In other words, p^ = L^d^ where the projection : I/-^([0, 1],M) A^ is characterised 
by 

/ dx= I d{x) dx, i = 1, . . . , A^, (4.1) 

for d G ^^"^([0, 1],M). Let ( , ) denote the duality pair between I/-^([0, 1],M) and its dual 
and let T/v be the space spanned by the characteristic functions Xi, • • • , Xn on the intervals 
/i, . . . , /at. Then (4.1) can alternatively be written as 

(L^rf, v) = (rf, v) for aU v G T^. (4.2) 

Thus, L^d is the Petrov-Galerkin projection of d (with respect to the spaces A^ and T/v). 

In Section 4.1, for a given map g : [0,1] [0,1], we constructed an approximate 
invariant density d G A^ by requiring that 

/ d{x) dx = d{x) dx, i = 1, . . . , N. (4.3) 

The next Lemma shows that this is equivalent to computing d G A^ such that 

L^PJ = J, 

where P is the Frobenius-Perron operator, i.e. is a fixed point of the discretised Frobenius- 
Perron operator 

Pn ■■= L'^P. (4.4) 
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Lemma 4.1. A function d G satisfies (4-3) if and only if it is a fixed point of the 
discrete Frohenius- Perron operator P^. 

Proof If d satisfies (4.3) then (J, v) = {Pd^ v) for all v G T/v since the characteristic 
functions on the intervals li are a basis of T/y. By the definition of the projection in 
(4.2), [U^Pd^v) = (Pd^v) for all v G T/y. Combining these two equalities, we arrive at 

(L^PJ, v) = (J, v) for aU v G Tn, 

and since d G and the projection is unique, it follows that P^d = d. 

To prove the converse, let d G A^ be a fixed point of P^, then {U^Pd^ v) = (rf, v) for 
all V G T/y. From the definition of L^, {U^Pd^v) = (Pd^v) for all v G T/y. Combining 
these two equalities, we obtain 

(PJ, ^) = (J, ^) for aU V eTN, 

which is equivalent to (4.3). □ 

Analogous to the proof of Lemma 8 in [11] we can show: 

Lemma 4.2. The discretised Frohenius- Perron operator P^ has a nonzero fixed point 

pn^n _ jn ^ An 

which satisfies 

d^(x) dx = 1. (4.5) 



We use the framework of [11] in order to prove convergence of our scheme in the case 
that g is piecewise C^[0, 1] and stretching (i.e. infa^^jo^i] Is'X^)! ^ The only difference 
to the setup in [11] is that we are dealing with a Petrov-Galerkin projection instead of a 
standard Galerkin projection. We are thus working with a fixed conjugate basis {^^}o of 

An 

A' = i = 0, . . . ,n, 

and choose the basis {a^jg of ^n+i such that 

{a\A^) = Sik, i,/c = 0, ...,n. 

Lemma 10 of [11] now yields stability of the projection L^: 

Lemma 4.3. There exist constants 71 (n) and 72 (^) such that for d G P^(0, 1) 

||L^rf||i<7i(n)|M||i 

and ifde BV{0, 1) then 

1 1 

\J Lid < ^2{n)\J d. (4.6) 
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As in [11], consistency of the projection is a standard result from approximation theory: 
Lemma 4.4. Let d G 1/^(0, 1), then 

hm ||d-L^rf||i = 0. 

If g is piecewise (7^ [0,1] and stretching, then the Lasota-Yorke theorem [33] imphes 
that for any function d G 1/^(0, 1) of bounded variation, Pd is of bounded variation as 
weU and there exist constants a > and /3 < 2/M (with M = 'mixe[o,i] \9\x)\ > 1) such 
that 

1 1 

\/ Pd<a\\d\\i + f3\/d. (4.7) 



Corollary 3 and 4 of [11] now yield convergence of our scheme: 

Theorem 4.5. Let g be piecewise (7^[0, 1] and stretching. If /3 < 1 and 72(^)/5 < 1 
then the sequence {d^)N of fixed points of the discretised Frohenius- Perron operator has a 
subsequence which converges in I/-^(0, 1) to an invariant density d of g. Moreover, 

\\d - rf^lli = 0{\\d - L^rflli) as N ^oo. 

Corollary 4.6. Under the conditions of Theorem 4-5 and for a given value of n, 

\\d -dlh = 0{h-+') = O (^^) asN^cx>. (4.8) 

Proof. The function L^d is the piecewise polynomial approximation to d assuming that 
the measures rrii on the intervals /j, i = 1, . . . , n are known exactly. This is precisely the 
situation that we considered in Theorem 3.1. Thus, 



\d - L^d||i = J \d{x) - dx 

N 

= J2 l^(^) ~ iL'}^d){x)\ dx 



1=1 

N 



J] / \D,{t) - PnAt)\ dt, 

i=l 



for appropriate values of ti and tr and where Di{t) and Pn,i{t) are obtained from the 
restriction of d{x) and L^d{x) respectively to the interval By (3.2), the final integral 
above is 0{h^^'^) and since a sum of = 0{h~^) such terms is required, then 

\\d-Lld\\^ = 0{h^+^). 

Finally, Theorem 4.5 tells us that 

\\d - d^lli = 0{\\d - L^dlli) = 0(/i"+i) = O 



Nn+l 

as required. □ 
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Figure 2: The map gi defined by (4.9) and it's density function di. 



Remark 4.7. Theorem 4.5 requires the map g to be stretching, i.e. infa;^[o,i] > 1 as 

this is a condition of the Lasota-Yorke Theorem [33]. However, by Theorem 3 of [33], this 
condition can be relaxed by requiring that g^ be stretching for some positive integer m 
and that infa^^[o,i] Is'X^)! ^ 0- Under these conditions. Theorem 4.5 also holds. 



4.3 Examples 

As an example of the preceding theory, we consider in detail the map ^'i : [0, 1] ^ [0, 1] 
defined by 

^, 0<x<V2-l 

9ii^)={ \ (4.9) 



^/2-l<x <l. 



2x 

The invariant density for this map (which is the map of [11]) is 

The function gi and the density di are shown in Fig. 2. 

We note that gi is piecewise C^[0, 1], but is not stretching as 5^1(1) = —1. However, 
g\ is stretching and infa;^[o,i] |5'i(^)| = 1 and so the conditions of Remark 4.7 are satisfied. 

We divided the interval [0, 1] into subintervals of width h = 1/N and took groups 
of n + 1 subintervals to construct a measure-preserving polynomial approximation to 
the density on these intervals. The method described in the previous section was used to 
compute the matrix representation A of the Frobenius-Perron operator and the eigenvector 
m of A associated with the eigenvalue nearest to +1 was computed using the inverse power 
method [5]. Since the true density function is known for this example, the one norm of the 
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10^ 



N 



10^ 



Figure 3: L^-errors in computed invariant densities for different values of n and for the 
map gi. 

difference between the exact function di{x) and the piecewise polynomial approximation 
d^{x) was determined. The program was written in Maple so that the accuracy of the 
calculations could be increased if required. The results are shown in Fig. 3. We note 
that the slope of the final segment of each of these lines is given for n = 0, 1, 2, 3 by 
— 1.035, —2.049, —3.049 and —3.917 respectively, which gives experimental verification of 
the theoretical result of Corollary 4.6 that the rate of convergence is 0{h^^^). 

To illustrate the power of this method, taking 16 subintervals with groups of 4 subin- 
tervals used to compute four cubic polynomial approximations to the density gives an 
error of 

\\di{x) - = 8.776219 x 10"^ 

Thus, by using higher order polynomials, very accurate results can be obtained using only 
a small number of subintervals. The standard Ulam method using only piecewise constant 
approximations to the density and the same number of intervals gives an error of 

\\di{x) - d%{x)\\i = 1.168727 x 10"^ 

For piecewise constant approximations, doubling approximately halves the error, as 
can be seen in Figure 3. Taking the error for = 128 and successively halving it, we 
obtain an approximation to the error for = 128 x 2'^ = 16, 384 of 1.307149 x 10"^/2'^ = 
1.02121 X 10~^, which is still slightly larger than the error above using pIq{x)\ 

We note that the results of [11] for this map appear to show convergence rates of 0(/i), 
0{h^) and 0{h^) for the piecewise constant, linear and quadratic cases. We get the same 
rates for the constant and linear cases, but we have 0{h^) for the piecewise quadratic 
case. 
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We have also applied our method to the map 



2x 1 

' 0<x<| 

92{x)={ (4.10) 

which is the map 6^4 of [11]. This map has invariant density given by 

d2[x) 



;i + x)2- 

Our results for this map are similar to those for the previous map and show the same 
rates of convergence. 
Finally, we consider 





1 












""2 






^3(^) 

which is the map S2 of [11]. In this case, the invariant density is 

d'^ix) = 12 I X 



which is quadratic. From the proof of Theorem 3.1, it is clear that the error term for the 
integrals involved in setting up the Frobenius-Perron operator depend on the (n + 1)*^ 
derivative of the density. Since the density is a quadratic function in this case, we would 
expect the error to be zero with a piecewise quadratic approximation and this is indeed 
the case. This is in contrast to the results of [11], whose results appear to show 0{h^) 
convergence with piecewise quadratic approximations for this map. We again achieve 
convergence rates of 0{h) and 0{h^) for piecewise constant and linear approximations. 

The interesting aspect of this example is that the function does not fulfil the condi- 
tions of Theorem 4.5, as it is neither C^[0, 1] (as g''^{x) = oc at the two points x = |±2~^/^) 
nor stretching (since the first and second derivatives at x = 1/2 are zero). Also no iterates 
of are stretching and so Remark 4.7 is no help either in this case. Thus, it seems that 
the method converges for a wider class of maps than those specified in Theorem 4.5 and 
Remark 4.7. In fact, recently much progress has been made in establishing the existence 
of invariant densities for general interval maps which are not expanding, see for example 
[1, 3, 4]. However, it remains to be explored whether the theory developed in these pa- 
pers yields the tools in order to prove Ulam's conjecture or the convergence of the scheme 
developed in this paper. 
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5 Computing the Lyapunov Exponent Using Integra- 
tion 

Having obtained a good approximation to the invariant density, we now want to compute 
an approximation to the Lyapunov exponent for the map g which is given by 



cr = / log \g\x)\d{x) dx. 
Jo 



Let d^ be the discontinuous piecewise polynomial approximation to the invariant density 
which is a nonzero fixed point of the discretised Frobenius-Perron operator given by 
(4.4). Then the approximation to the Lyapunov exponent that we can actually compute 
(ignoring numerical integration errors) is given by 



(Jn= I log|^'(x)|rf^(x) dx. 
Jo 



A simple error analysis gives that 



log \g\x)\{d{x) — d^{x)) dx 



< [ \log\g\x)\{d{x) - d^'^ix))] dx 
Jo 

< \\logW\ \\^\\d-d%\\, 
= 0(/i"+^), 

using Holder's inequality and (4.8), assuming that logl^''! is bounded. Thus, it would 
appear that the error in the Lyapunov exponent is determined by the error in the invariant 
density. However, we now show that better results than this can often be obtained. 

Theorem 5.1. Assume that g is piecewise C^[0, 1] and stretching with /3 < 1 (cf. (4-7)) 
and 72(77.) ;5 < 1 (cf. (4-6)), and that \g'\ G (7°[0, 1]. If the Frobenius-Perron operator P 
has a unique invariant density d G (7^^^[0, 1] (n odd) or d G C^^^[0^ 1] (n even) then 



n+2\ 



as N = 1/h 00. 



Proof. We consider two different approaches to proving this result. 
For the first proof, we note from Corollary 4.6 that 

and so we write 

d{x) =rf^(x) + /i^+^5(x). (5.1) 
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where 5 G L^{0, 1). Integrating by parts then gives 

»i 

log \g' {x)\(d{x) — d^{x)) dx 
\og\g'{x)\ I d{x') - d'l^^x') dx' 



.^n+l j i_lQg|^/(2.)| r 8{x') dx' dx. 

Jo d^ Jo 



Clearly, the boundary term evaluated at x = is zero and since 



d{x) dx 



dj\^(^x^ dx — 1, 



by (4.5), then the boundary term evaluated at x = 1 is also zero. Using the same methods 
as in the proof of Theorem 3.3, this integration by parts results in the error increasing by 
one power of h and so the error is 0(/i^^^) as claimed. 

For the second proof, we let = U^d be the projection of the true invariant density 
d onto the approximation space A^, where the projection is defined by (4.1). Then 

\og\g' {x)\{d{x) -p%{x)) dx 



\^-^n\ = 



< 



log |^'(x)|(p^(x) - rf^(x)) dx 
log\g\x)\{d{x) -p%{x)) dx 
log \g' {x)\{p%{x) - rf^(x)) dx 



The first integral is of the type that we considered in Section 3.2. The stated conditions 
on g ensure that log Is'X^)! satisfies the conditions of Theorem 3.4 and so this Theorem 



gives 



\og\g' {x)\{d{x) -p%{x)) dx = 



0(/i^+3), nodd 



0{h 



n+2) 



n even 



For the second integral, we first note that 

Pn ~ ~ L^{d — d^) 

using (5.1), where 



ft Ojy^ 



s:n _ jn ^ 



Integrating the second integral by parts then gives 

/ \og\g' {x)\{jf'j^{x) - d:!l;{x)) dx = \og\g'{x)\ p^^^x') - d^^^^x') dx' 
Jo I Jo Jo 

/" ^\og\g'{x)\ ! 5l{x')dx'dx. 
Jo d.x Jq 
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Figure 4: Absolute errors in the computed values of the Lyapunov exponent for the map 
gi{x). 

Again, the boundary term at x = is zero and since 

/ p^{x) dx = d^{x) dx = 1, 
Jo Jo 

by (4.5), then the boundary term evaluated at x = 1 is also zero. Using the same approach 
now as for the first proof above, the error has increased by one power of h and so this 
error is 0(/i^+^). 

Combining these two integrals, asymptotically the largest error is from the second 
integral and so the error in the Lyapunov exponent is 0(/i^^^) as claimed. □ 

5.1 Example 

We consider again the map gi defined by (4.9) that we considered previously in Section 
4.3. We note that for this example, |5'i(x)| is continuous at the point x = ^/2 — 1, which 
is one of the conditions of Theorem 5.1. 

We used the method described in Section 4 to obtain the piecewise polynomial ap- 
proximation rf^(x) to the density d{x) and we then evaluated 

(Jn= I log|^;(x)|rf^(x) dx. 
Jo 

The calculations were done in Maple and so the integration was performed accurately 
by a Maple routine rather than using Gaussian quadrature. 

Since we know the true density function for this map, we can also accurately compute 
the true value of the Lyapunov exponent a, which is found numerically to have the value 
log 2. The results of these computations are shown in Fig. 4. 

The slope of the final segment of the lines for n = and n = 2 are given by —2.031 and 
—4.001 respectively, and so, for even values of n, the results agree well with the prediction 
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of Theorem 5.1 that the rate of convergence is 0(/i^+^). However, for n = 1, the rate 
of convergence initiaUy is very similar to that for n = 2, which is 0{h^)^ or equivalently 
0(/i^+^), and this is higher than that predicted by Theorem 5.1. However, we note that 
the slope of the final segment of the line for n = 1 is —3.105, which is much closer to 
the predicted rate of convergence of 0{h^). Similarly, for n = 3, the rate of convergence 
estimate initially oscillates, but seems to be higher than the prediction of 0{h^). However, 
considering again the slope of the final segment of the line, we obtain a value of —5.091, 
which is close to the predicted rate of 0{h^). 

These results can be understood from the second proof of Theorem 5.1. In that proof, 
the error was broken up into two terms. For n even, the two separate terms were both 
0(/i^+^). However, for n odd, the terms were of diflFerent order so that approximately 

In this case, the asymptotic rate of convergence as /i ^ is clearly However, the 

rate of convergence that is observed for finite values of h also depends on the magnitude 
of the two constants Ci and C2. If Ci 3> C2, then for relatively "large" values of /i, the 
dominant term in the error will be the first one and so the rate of convergence will appear 
to be 0{h^^^). However, as h decreases, eventually the higher power of h will ensure that 
the first term becomes smaller than the second, and then the true rate of convergence of 
0(/i^^^) will be observed. Clearly, this is what is happening in our example. 

We note that even though the n = 1 and n = 2 cases have diflFerent asymptotic rates of 
convergence, the actual magnitude of the errors shown in Figure 4 are similar in these two 
cases. We have noted that for odd n, the first integral which is 0{h^^^) seems to dominate 
the errors for moderate values of /i, and this implies that the error initially decreases at 
a similar rate to that for the next even value of n. Thus, even though theoretically the 
results for n = 2 should be better than those forn = 1, and would be better for sufficiently 
small /i, for moderate values of /i, the results for these two values are comparable. 

Similar results are obtained for the map g2{x) given by (4.10). 

6 Future Directions 
6.1 Stochastic perturbations 

By construction, our approach requires the invariant density to be smooth. For determin- 
istic maps, this is not the generic situation. However, in applications one is often faced 
with a system which is additionally perturbed by (small) random inffuences. In these 
cases, instead of the deterministic system one considers the dynamics 

= g{x) 

where ^ is chosen from a given probability distribution /x. Suppose that /x is absolutely 
continuous with density h. Then the associated Frobenius-Perron operator on is given 

by 

Pf{x) = J h{g{y) - x)f{y) dy, f e L\ (6.1) 
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For common distributions (like a normal distribution or the uniform distribution sup- 
ported on a small ball around 0) the associated Markov chain possesses finitely many 
invariant measures which, according to (6.1), inherit the smoothness of the distribution 



In this setting, the approach proposed in this paper yields a highly accurate method 
for the approximation of the invariant distribution. For very small perturbations and 
in the case of a non-smooth invariant density, however, one will still need many modes 
for this. It would be interesting to investigate how the approximation error behaves in 
dependence on the magnitude of the perturbation. 

6.2 Higher dimensions 

The methods described above for one-dimensional maps can easily be generalised to higher 
dimensions. We briefiy consider only the case of a bilinear approximation to the density 
given the total measure on four neighbouring squares in two dimensions. We assume that 
a change of variables has been performed so that the region of interest is [—1,1]^. We 
take t and r as the two independent variables and we want to approximate the density 
function D{t^T) over this region. 

We first define to = "R) = — 1, = "T"! = and t2 = "^2 = 1 and assume that we know 
the four values 



We then want to construct a bilinear approximation to the density which preserves the 
total measure on each of the four subregions. We write the polynomial approximation as 

2 



This polynomial basis can then be used to approximate the invariant density, again giving 
a discontinuous approximation over the whole region. A similar error analysis as in the 
one-dimensional case can also be performed. 



[10]. 




k,i=i 
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